Observed log-chlorophyll at representative station for the St. Lucie Estuary

library(tidyverse)
library(lubridate)
library(mgcv)  
library(plotly)
library(WRTDStidal)
library(gridExtra)
source('R/funcs.R')

# format the data to model
data(sl_dat)
sl_mod <- sl_dat %>%
  rename(date = Date) %>% 
  mutate(
    doy = yday(date), 
    dec_time = decimal_date(date), 
    yr = year(date),
    mo = month(date, label = T)
  ) %>% 
  mutate(
    flo = sal, 
    lnchl = log(1 + chl)
    ) %>% 
  select(-sal)

# plot, all
p <- ggplot(sl_mod, aes(x = date, y = lnchl)) + 
  geom_line() +
  theme_bw() 
ggplotly(p)
# boxplot, by years
p <- ggplot(sl_mod, aes(x = yr, y = lnchl)) + 
  geom_boxplot() + 
  theme_bw()
ggplotly(p)
# boxplot, by month
p <- ggplot(sl_mod, aes(x = mo, y = lnchl)) + 
  geom_boxplot() + 
  theme_bw()
ggplotly(p)

Some simple GAMs to explore annual, seasonal trends.

# annual only
mod1 <- gam(lnchl ~ s(dec_time, bs = 'tp'),
  data = sl_mod, 
  na.action = na.exclude
  )

# seasonal only
mod2 <- gam(lnchl ~ s(doy, bs = 'cc'), 
  knots = list(doy = c(1, 366)),
  data = sl_mod, 
  na.action = na.exclude
  )
 
# annual and seasonal, additive
mod3 <- gam(lnchl ~ s(dec_time, bs = 'tp') + 
  s(doy, bs = 'cc'), 
  knots = list(doy = c(1, 366)),
  data = sl_mod, 
  na.action = na.exclude
  )

# annual and seasonal, interaction
mod4 <- gam(lnchl ~ te(dec_time, doy, bs = c('tp', 'cc')),
  knots = list(doy = c(1, 366)),
  data = sl_mod, 
  na.action = na.exclude
  )

Summary stats of annual, seasonal models:

mods <- list(
  mod1 = mod1,
  mod2 = mod2, 
  mod3 = mod3,
  mod4 = mod4
  )

# smoother stats of GAMs
map(mods, ~ summary(.x)$s.table %>% data.frame %>% rownames_to_column('smoother')) %>% 
  enframe %>% 
  unnest %>% 
  kable
name smoother edf Ref.df F p.value
mod1 s(dec_time) 1.646162 2.048370 9.703867 6.78e-05
mod2 s(doy) 4.098673 8.000000 6.815467 0.00e+00
mod3 s(dec_time) 1.943078 2.429282 8.906536 5.67e-05
mod3 s(doy) 4.114540 8.000000 7.052094 0.00e+00
mod4 te(dec_time,doy) 9.994763 12.933682 5.562618 0.00e+00 
# summary stats of GAMs
map(mods, ~ data.frame(
    AIC = AIC(.x), 
    R2 = summary(.x)$r.sq)) %>% 
  enframe %>% 
  unnest %>% 
  kable
name AIC R2
mod1 683.2360 0.0500261
mod2 653.5304 0.1302894
mod3 634.5680 0.1787925
mod4 643.6646 0.1668405 
pred_dat <- sl_mod

# predictions
sl_res <- pred_dat %>% 
  select(date, flo, dec_time, doy, yr) %>% 
  mutate(
    mod1 = predict(mod1, newdata = pred_dat), 
    mod2 = predict(mod2, newdata = pred_dat), 
    mod3 = predict(mod3, newdata = pred_dat), 
    mod4 = predict(mod4, newdata = pred_dat)
  ) %>% 
  tidyr::gather('mods', 'pred', -date, -doy, -dec_time, -flo, -yr)

# plot
p <- ggplot(sl_res, aes(x = date)) + 
  geom_point(data = sl_mod, aes(y = lnchl), size = 0.5) + 
  geom_line(aes(y = pred, colour = mods)) + 
  theme_bw() + 
  theme(
    legend.position = 'top', 
    legend.title = element_blank()
    )
ggplotly(p)
# plot
p <- ggplot(sl_res, aes(x = doy, group = factor(yr), colour = yr)) + 
  geom_line(aes(y = pred)) + 
  theme_bw() + 
  theme(
    legend.position = 'top', 
    legend.title = element_blank()
    ) + 
  facet_wrap(~ mods, ncol = 2)
ggplotly(p)

Adding flow data to the model:

# model, all terms additive
mod5 <- gam(lnchl ~ s(dec_time, bs = 'tp') + s(doy, bs = 'cc') + s(flo, bs = 'tp'),
  knots = list(doy = c(1, 366)),
  data = sl_mod,
  na.action = na.exclude
  )

# model, flo additive, interaction between dec_time, doy
mod6 <- gam(lnchl ~ te(dec_time, doy, bs = c('tp', 'cc')) + s(flo, bs = 'tp'),
  knots = list(doy = c(1, 366)),
  data = sl_mod,
  na.action = na.exclude
  )

# model, doy additive, interaction between dec_time and flow
mod7 <- gam(lnchl ~ te(dec_time, flo, bs = c('tp', 'tp')) + s(doy, bs = 'cc'),
  knots = list(doy = c(1, 366)),
  data = sl_mod,
  na.action = na.exclude
  )

# model, dec_time additive, interaction between flo and doy
mod8 <- gam(lnchl ~ te(flo, doy, bs = c('tp', 'cc')) + s(dec_time, bs = 'tp'),
  knots = list(doy = c(1, 366)),
  data = sl_mod,
  na.action = na.exclude
  )

# model, interaction between flow and doy, interaction between flo and dec_time
mod9 <- gam(lnchl ~ te(flo, doy, bs = c('tp', 'cc')) + te(flo, dec_time, bs = c('tp', 'tp')),
  knots = list(doy = c(1, 366)),
  data = sl_mod,
  na.action = na.exclude
  )

# model, all terms interaction
mod10 <- gam(lnchl ~ te(dec_time, doy, flo, bs = c('tp', 'cc', 'tp')),
  knots = list(doy = c(1, 366)),
  data = sl_mod,
  na.action = na.exclude
  )

Summary stats of annual, seasonal, flow models:

mods2 <- list(
  mod5 = mod5,
  mod6 = mod6, 
  mod7 = mod7,
  mod8 = mod8,
  mod9 = mod9, 
  mod10 = mod10
  )

# smoother stats of GAMs
map(mods2, ~ summary(.x)$s.table %>% data.frame %>% rownames_to_column('smoother')) %>% 
  enframe %>% 
  unnest %>% 
  kable
name smoother edf Ref.df F p.value
mod5 s(dec_time) 2.024919 2.531950 7.173692 0.0003743
mod5 s(doy) 3.670840 8.000000 6.512722 0.0000000
mod5 s(flo) 2.821335 3.513840 4.559861 0.0025490
mod6 te(dec_time,doy) 7.580584 9.991352 6.465608 0.0000000
mod6 s(flo) 2.843758 3.538719 4.837179 0.0016166
mod7 te(dec_time,flo) 9.293442 11.920524 4.107297 0.0000058
mod7 s(doy) 3.634232 8.000000 6.279981 0.0000000
mod8 te(flo,doy) 7.669868 10.031915 7.137935 0.0000000
mod8 s(dec_time) 2.090841 2.610887 7.828485 0.0001493
mod9 te(flo,doy) 8.154234 10.577257 6.447297 0.0000000
mod9 te(flo,dec_time) 6.190206 20.000000 1.310943 0.0000112
mod10 te(dec_time,doy,flo) 27.336212 37.654333 3.674406 0.0000000 
# summary stats of GAMs
map(mods2, ~ data.frame(
    AIC = AIC(.x), 
    R2 = summary(.x)$r.sq)) %>% 
  enframe %>% 
  unnest %>% 
  kable
name AIC R2
mod5 527.8769 0.2291574
mod6 534.5420 0.2171463
mod7 522.0281 0.2536808
mod8 525.3183 0.2383925
mod9 522.0352 0.2568790
mod10 510.9187 0.3099731 
# predictions
sl_res2 <- pred_dat %>% 
  select(date, flo, dec_time, doy, yr) %>% 
  mutate(
    mod5 = predict(mod5, newdata = pred_dat), 
    mod6 = predict(mod6, newdata = pred_dat), 
    mod7 = predict(mod7, newdata = pred_dat), 
    mod8 = predict(mod8, newdata = pred_dat),
    mod9 = predict(mod9, newdata = pred_dat),
    mod10 = predict(mod10, newdata = pred_dat)
  ) %>% 
  tidyr::gather('mods', 'pred', -date, -doy, -dec_time, -flo, -yr)

# plot
p <- ggplot(sl_res2, aes(x = date)) + 
  geom_point(data = sl_mod, aes(y = lnchl), size = 0.5) + 
  geom_line(aes(y = pred, colour = mods)) + 
  theme_bw() + 
  theme(
    legend.position = 'top', 
    legend.title = element_blank()
    )
ggplotly(p)
ptheme <- theme(
  axis.title.x = element_blank(), 
  axis.title.y = element_blank()
)
cols <- 'Spectral'
p5 <- dynagam(mod5, pred_dat, ncol = 1, col_vec = cols) + 
  ptheme + 
  theme(legend.position = 'none') +
  ggtitle('mod5')
p6 <- dynagam(mod6, pred_dat, ncol = 1, col_vec = cols) + 
  ptheme + 
  theme(legend.position = 'none') + 
  ggtitle('mod6')
p7 <- dynagam(mod7, pred_dat, ncol = 1, col_vec = cols) + 
  ptheme + 
  theme(legend.position = 'none') + 
  ggtitle('mod7')
p8 <- dynagam(mod8, pred_dat, ncol = 1, col_vec = cols) + 
  ptheme + 
  theme(legend.position = 'none') + 
  ggtitle('mod8')
p9 <- dynagam(mod9, pred_dat, ncol = 1, col_vec = cols) + 
  ptheme + 
  theme(legend.position = 'none') + 
  ggtitle('mod9')
p10 <- dynagam(mod10, pred_dat, ncol = 1, col_vec = cols) + 
  ptheme + 
  ggtitle('mod10')
pleg <- g_legend(p10)
p10 <- p10 + 
  theme(legend.position = 'none')

grid.arrange(
  pleg, 
  arrangeGrob(p5, p6, p7, p8, p9, p10, ncol = 6, bottom = 'lnQ', left = 'lnchl'), 
  ncol = 1, 
  heights = c(0.1, 1)
)

Backwards model selection, see here:

mod <- gam(lnchl ~ s(dec_time, bs = 'tp') + 
              s(doy, bs = 'cc') + 
              s(flo, bs = 'tp') +
              te(flo, doy, bs = c('tp', 'cc')) + 
              te(flo, dec_time, bs = c('tp', 'tp')) +
              te(dec_time, doy, bs = c('tp', 'cc')) +
              te(dec_time, doy, flo, bs = c('tp', 'cc', 'tp')),
            knots = list(doy = c(1, 366)),
            data = sl_mod,
            na.action = na.exclude
)
summary(mod)
AIC(mod)

mod <- gam(lnchl ~ s(dec_time, bs = 'tp') + 
              s(doy, bs = 'cc') + 
              s(flo, bs = 'tp') +
              te(flo, doy, bs = c('tp', 'cc')) + 
              te(flo, dec_time, bs = c('tp', 'tp')) +
              te(dec_time, doy, flo, bs = c('tp', 'cc', 'tp')),
            knots = list(doy = c(1, 366)),
            data = sl_mod,
            na.action = na.exclude
)
summary(mod)
AIC(mod)

mod <- gam(lnchl ~ s(dec_time, bs = 'tp') + 
              s(doy, bs = 'cc') + 
              te(flo, doy, bs = c('tp', 'cc')) + 
              te(flo, dec_time, bs = c('tp', 'tp')) +
              te(dec_time, doy, flo, bs = c('tp', 'cc', 'tp')),
            knots = list(doy = c(1, 366)),
            data = sl_mod,
            na.action = na.exclude
)
summary(mod)
AIC(mod)

mod <- gam(lnchl ~ s(dec_time, bs = 'tp') + 
              s(doy, bs = 'cc') + 
              te(flo, dec_time, bs = c('tp', 'tp')) +
              te(dec_time, doy, flo, bs = c('tp', 'cc', 'tp')),
            knots = list(doy = c(1, 366)),
            data = sl_mod,
            na.action = na.exclude
)
summary(mod)
AIC(mod)

mod <- gam(lnchl ~ s(dec_time, bs = 'tp') + 
              s(doy, bs = 'cc') + 
              te(flo, doy, bs = c('tp', 'cc')) + 
              te(dec_time, doy, flo, bs = c('tp', 'cc', 'tp')),
            knots = list(doy = c(1, 366)),
            data = sl_mod,
            na.action = na.exclude
)
summary(mod)
AIC(mod)

mod <- gam(lnchl ~ s(dec_time, bs = 'tp') + 
              s(doy, bs = 'cc') + 
              te(dec_time, doy, flo, bs = c('tp', 'cc', 'tp')),
            knots = list(doy = c(1, 366)),
            data = sl_mod,
            na.action = na.exclude
)
summary(mod)
AIC(mod)

mod <- gam(lnchl ~ s(doy, bs = 'cc') + 
              te(dec_time, doy, flo, bs = c('tp', 'cc', 'tp')),
            knots = list(doy = c(1, 366)),
            data = sl_mod,
            na.action = na.exclude
)
summary(mod)
AIC(mod)